Requirements

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(readr)
library(tf)
## 
## Attaching package: 'tf'
## 
## The following objects are masked from 'package:stats':
## 
##     sd, var
library(tidyfun)
library(ggplot2)
library(refund)
library(mgcViz)
## Loading required package: mgcv
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## 
## The following object is masked from 'package:dplyr':
## 
##     collapse
## 
## This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
## Loading required package: qgam
## Registered S3 method overwritten by 'mgcViz':
##   method from   
##   +.gg   ggplot2
## 
## Attaching package: 'mgcViz'
## 
## The following objects are masked from 'package:stats':
## 
##     qqline, qqnorm, qqplot
library(rgl)
library(plotly)
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout
library(fdaoutlier)
library(purrr)
library(patchwork)
library(table1)
## 
## Attaching package: 'table1'
## 
## The following objects are masked from 'package:base':
## 
##     units, units<-
sessionInfo()
## R version 4.5.1 (2025-06-13)
## Platform: x86_64-pc-linux-gnu
## Running under: Linux Mint 22.1
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Berlin
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] table1_1.4.3     patchwork_1.3.0  fdaoutlier_0.2.1 plotly_4.10.4   
##  [5] rgl_1.3.18       mgcViz_0.2.0     qgam_2.0.0       mgcv_1.9-3      
##  [9] nlme_3.1-168     refund_0.1-37    tidyfun_0.0.98   tf_0.3.5        
## [13] lubridate_1.9.4  forcats_1.0.0    stringr_1.5.1    dplyr_1.1.4     
## [17] purrr_1.0.4      readr_2.1.5      tidyr_1.3.1      tibble_3.3.0    
## [21] ggplot2_3.5.2    tidyverse_2.0.0 
## 
## loaded via a namespace (and not attached):
##   [1] Rdpack_2.6.4       bitops_1.0-9       gridExtra_2.3     
##   [4] rlang_1.1.6        magrittr_2.0.3     matrixStats_1.5.0 
##   [7] compiler_4.5.1     vctrs_0.6.5        pkgconfig_2.0.3   
##  [10] fastmap_1.2.0      backports_1.5.0    magic_1.6-1       
##  [13] promises_1.3.2     deSolve_1.40       rmarkdown_2.29    
##  [16] tzdb_0.5.0         pracma_2.4.4       nloptr_2.2.1      
##  [19] xfun_0.52          cachem_1.1.0       jsonlite_2.0.0    
##  [22] later_1.4.2        parallel_4.5.1     cluster_2.1.8.1   
##  [25] R6_2.6.1           bslib_0.9.0        stringi_1.8.7     
##  [28] RColorBrewer_1.1-3 GGally_2.3.0       extrafontdb_1.0   
##  [31] boot_1.3-31        rainbow_3.8        jquerylib_0.1.4   
##  [34] Rcpp_1.0.14        iterators_1.0.14   knitr_1.50        
##  [37] zoo_1.8-14         base64enc_0.1-3    hdrcde_3.4        
##  [40] extrafont_0.19     httpuv_1.6.16      Matrix_1.7-3      
##  [43] splines_4.5.1      timechange_0.3.0   tidyselect_1.2.1  
##  [46] rstudioapi_0.17.1  abind_1.4-8        yaml_2.3.10       
##  [49] viridis_0.6.5      doParallel_1.0.17  codetools_0.2-20  
##  [52] miniUI_0.1.1.1     RLRsim_3.1-8       lattice_0.22-5    
##  [55] plyr_1.8.9         shiny_1.10.0       ks_1.14.3         
##  [58] withr_3.0.2        S7_0.2.0           pbs_1.1           
##  [61] evaluate_1.0.3     ggstats_0.9.0      grpreg_3.5.0      
##  [64] mclust_6.1.1       pillar_1.10.2      fda_6.2.0         
##  [67] KernSmooth_2.23-26 checkmate_2.3.2    foreach_1.5.2     
##  [70] reformulas_0.4.1   pcaPP_2.0-5        generics_0.1.3    
##  [73] RCurl_1.98-1.17    hms_1.1.3          scales_1.4.0.9000 
##  [76] minqa_1.2.8        xtable_1.8-4       gamm4_0.2-7       
##  [79] glue_1.8.0         lazyeval_0.2.2     tools_4.5.1       
##  [82] fds_1.8            data.table_1.17.0  lme4_1.1-37       
##  [85] mvtnorm_1.3-3      grid_4.5.1         Rttf2pt1_1.3.12   
##  [88] rbibutils_2.3      colorspace_2.1-1   Formula_1.2-5     
##  [91] cli_3.6.5          viridisLite_0.4.2  gtable_0.3.6      
##  [94] sass_0.4.10        digest_0.6.37      htmlwidgets_1.6.4 
##  [97] farver_2.1.2       htmltools_0.5.8.1  lifecycle_1.0.4   
## [100] httr_1.4.7         mime_0.13          MASS_7.3-65
strain_df<-readRDS("papillary muscle_data.rds")
strain_df$id<-factor(strain_df$id)
scalar_covariates<-readRDS("scalar_covariates.rds")
scalar_covariates$id<-factor(scalar_covariates$id)
# put irregular functions x that vary in length/start/end on a new_domain
# on a regular grid of length n_arg
tf_rectify <- function(x, new_domain = c(0, 1), n_arg = 30, ...) {
  x_val <- tf_evaluations(x)
  x_arg <- tf_arg(x) |> ensure_list()
  # create new grids with same proportional spacings as old grids on the new domain:
  new_arg_irreg <- purrr::map(x_arg, 
                        function(x) (((x - min(x))/diff(range(x))) * # scale to [0,1]
                                        diff(new_domain)) +  # scale to [0, new_length]
                                       new_domain[1]) # push to new domain
  new_arg_grid <- seq(new_domain[1], new_domain[2], length = n_arg)
  # put data on new common domain and re-eval on new grid
  tfd(x_val, new_arg_irreg, ...) |> tfd(arg = new_arg_grid, ...) 
}

Descriptives / Data Viz / Data Prep

# get regular evenly spaced time points for the curves
strain_df_clean <- strain_df %>%
  mutate(
    strain_str = ifelse(strain_curve %in% c("GLS4", "GLS2", "GLS3"), "GLS",
                        ifelse(strain_curve %in% c("LA4", "LA2"), "LA", strain_curve)
    ),
    strain_reg = tfd(strain_val[, seq(0, 1, length.out = 60), interpolate = TRUE]),
    id = factor(id)
  )

# show complete dataset before aggregation:
ggplot(strain_df_clean) +
  geom_spaghetti(aes(y = strain_reg, color = strain_str, linetype = cycle_number)) +
  facet_wrap(~id) + labs(title = "Full Dataset")

# global mean strain curves
strain_df_clean |>
  group_by(strain_str) %>%
  summarize(
    strain_mean = mean(strain_reg),
    strain_se = sd(strain_reg) / sqrt(nlevels(id))
  ) %>%
  filter(strain_str %in% c("APM", "PPM", "TOR", "LA", "GLS")) %>%
  ggplot() +
  geom_spaghetti(aes(y = strain_mean, color = strain_str)) +
  geom_errorband(aes(
    ymin = strain_mean - 2 * strain_se,
    ymax = strain_mean + 2 * strain_se,
    fill = strain_str
  )) + labs(title = "Global Mean Strain Curves")

# average curves per patient: gls=mean(gls2,gls3, gls4), LA=mean (LA2, LA4)
strain_df_patientwise <- strain_df_clean %>%
  group_by(id, strain_str) %>%
  summarise(
    strain_m = mean(strain_reg)
  )
## `summarise()` has grouped output by 'id'. You can override using the `.groups`
## argument.
ggplot(strain_df_patientwise) +
  geom_spaghetti(aes(y = strain_m, color = strain_str)) +
  facet_wrap(~id) + labs(title = "Patient-wise Mean Strain Curves")

# putting in right format for analysis and  adding the scalar covariates of interest
strain_df_patientwise_long <- strain_df_patientwise %>%
  dplyr::select(c("id", "strain_str", "strain_m")) %>%
  pivot_wider(names_from = strain_str, values_from = strain_m)

strain_df_patientwise_long <- strain_df_patientwise_long[-43, ]

strain_df_patientwise_long <- left_join(strain_df_patientwise_long, scalar_covariates, by = "id") #
strain_analysisdata <- strain_df_patientwise_long

Table 1

strain_analysisdata_summaries <- strain_analysisdata %>% 
    mutate(
      tor_max_time = TOR|>tf_where(value == max(value), "first"),
         tor_max = TOR|>tf_fmax(),
         apm_min_time = APM|>tf_where(value == min(value), "first"),
         apm_min = APM|>tf_fmin(),
         ppm_min_time=PPM|>tf_where(value == min(value), "first"),
         ppm_min=PPM|>tf_fmin(),
         gls_min_time=GLS|>tf_where(value == min(value), "first"),
         gls_min=GLS |> tf_fmin())

table1(~ 
         Age + factor(Gender) + factor(HTN) + factor(DIABETES) + weight + height + 
         heart_rate + tor_max + tor_max_time + gls_min + gls_min_time + 
         apm_min + apm_min_time + ppm_min + ppm_min_time,
       data = strain_analysisdata_summaries)
Overall
(N=75)
Age
Mean (SD) 37.4 (14.9)
Median [Min, Max] 36.2 [15.6, 78.1]
factor(Gender)
0 16 (21.3%)
1 59 (78.7%)
factor(HTN)
0 63 (84.0%)
1 12 (16.0%)
factor(DIABETES)
0 68 (90.7%)
1 7 (9.3%)
weight
Mean (SD) 75.0 (13.3)
Median [Min, Max] 76.5 [47.0, 115]
Missing 5 (6.7%)
height
Mean (SD) 173 (8.21)
Median [Min, Max] 175 [151, 191]
Missing 5 (6.7%)
heart_rate
Mean (SD) 75.3 (13.7)
Median [Min, Max] 72.5 [43.0, 128]
Missing 1 (1.3%)
tor_max
Mean (SD) 25.2 (10.6)
Median [Min, Max] 24.0 [9.57, 54.1]
tor_max_time
Mean (SD) 0.412 (0.0649)
Median [Min, Max] 0.407 [0.271, 0.593]
gls_min
Mean (SD) -17.3 (2.20)
Median [Min, Max] -17.2 [-22.8, -12.2]
Missing 3 (4.0%)
gls_min_time
Mean (SD) 0.429 (0.0557)
Median [Min, Max] 0.441 [0.288, 0.593]
Missing 3 (4.0%)
apm_min
Mean (SD) -17.1 (3.80)
Median [Min, Max] -17.0 [-29.7, -7.22]
apm_min_time
Mean (SD) 0.509 (0.0765)
Median [Min, Max] 0.492 [0.322, 0.712]
ppm_min
Mean (SD) -17.8 (3.73)
Median [Min, Max] -17.9 [-25.4, -8.63]
ppm_min_time
Mean (SD) 0.509 (0.0791)
Median [Min, Max] 0.508 [0.322, 0.695] 
m3_data <- strain_analysisdata %>%
  select(id, APM, PPM, TOR, Age, Gender, heart_rate)

# plotting the apm, ppm and torsion curves
f_summary <- m3_data %>%
  dplyr::select("id", "APM", "PPM", "TOR") %>%
  pivot_longer(!id, names_to = "Curve", values_to = "strain") %>%
  group_by(Curve) %>%
  summarise(
    strain_mean = mean(strain),
    strain_sd = sd(strain)
  ) %>%
  mutate(
    maximum = strain_mean |> tf_where(value == max(value), "first"),
    minimum = strain_mean |> tf_where(value == min(value), "first")
  )


f_summary %>%
  ggplot() +
  geom_spaghetti(aes(y = strain_mean, color = Curve)) +
  geom_errorband(aes(
    ymax = strain_mean + 2 * strain_sd / sqrt(nrow(m3_data)),
    ymin = strain_mean - 2 * strain_sd / sqrt(nrow(m3_data)),
    fill = Curve
  )) +
  geom_vline(xintercept = 0.41, color = "blue", linetype = 2) + # time of maximum torsion
  geom_vline(xintercept = 0.50, linetype = 2) + # time of minimum PM (ie.maximum contraction)
  labs(
    x = "Cardiac cycle proportion",
    y = "",
    caption = "mean+/-2 standard errors",
    legend = "Curve"
  ) +
  theme_bw()

Peak Timing Difference Analysis

# comparing time of maximum PM contraction to time of maximum torsion occurrence
APM_min<-m3_data$APM |> tf_where(value == min(value), "first")
PPM_min<-   m3_data$PPM |> tf_where(value == min(value), "first")
Tor_max<-m3_data$TOR |> tf_where(value == max(value), "first")

mean(c(APM_min,PPM_min))
## [1] 0.509
mean(APM_min)
## [1] 0.509
sd(APM_min)
## [1] 0.0765
mean(PPM_min)
## [1] 0.509
sd(PPM_min)
## [1] 0.0791
mean(Tor_max)
## [1] 0.412
sd(Tor_max)
## [1] 0.0649
(res1 <- t.test(APM_min, Tor_max, paired = TRUE))
## 
##  Paired t-test
## 
## data:  APM_min and Tor_max
## t = 11, df = 74, p-value <2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.0798 0.1146
## sample estimates:
## mean difference 
##          0.0972
(res2 <- t.test(PPM_min, Tor_max, paired = TRUE))
## 
##  Paired t-test
## 
## data:  PPM_min and Tor_max
## t = 14, df = 74, p-value <2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.0828 0.1111
## sample estimates:
## mean difference 
##          0.0969
(res3 <- t.test(PPM_min, APM_min, paired = TRUE))
## 
##  Paired t-test
## 
## data:  PPM_min and APM_min
## t = -0.02, df = 74, p-value = 1
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.0183  0.0179
## sample estimates:
## mean difference 
##       -0.000226
# min apm/ppm ( maximum of papillary muscles contraction) occurs significantly  later than maximum torsion 

# association across patients:
layout(t(1:2))
plot(APM_min, Tor_max, asp = 1, col = rgb(0,0,0,.3)); abline(c(0,1), col = "gold") 
plot(PPM_min, Tor_max, asp = 1, col = rgb(0,0,0,.3)); abline(c(0,1), col = "gold")

layout(t(1:2))
plot(APM_min, Tor_max - APM_min, asp = 1, col = rgb(0,0,0,.3)); abline(c(0,0), col = "gold") 
plot(PPM_min, Tor_max - PPM_min, asp = 1, col = rgb(0,0,0,.3)); abline(c(0,0), col = "gold")

# patients with later APM_min (PPM_min as well to lesser extent) tend to have larger delay between Tor_max and APM_min
# selecting only diastolic curve value, diastole begins at the maximum of torsion
m4_data <- strain_analysisdata %>%
  select(id, TOR, APM, PPM, Gender, Age, heart_rate) %>%
  mutate(tor_max = TOR |> tf_where(value == max(value), "first")) %>%
  mutate(
    tor_diastole = TOR |> tf_zoom(tor_max, 1),
    apm_diastole = APM |> tf_zoom(tor_max, 1),
    ppm_diastole = PPM |> tf_zoom(tor_max, 1)
  )
## Warning: There were 222 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `tor_diastole = tf_zoom(TOR, tor_max, 1)`.
## Caused by warning:
## ! Combining incompatible <tfd_reg> with <tfd_reg> by casting to <tfd_irreg>.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 221 remaining warnings.
m4_data$heart_rate[is.na(m4_data$heart_rate)] <- 
  mean(m4_data$heart_rate, na.rm = TRUE) # imputing the mean for one missing heart rate observation

# obtaining regular and same domain over diastole torsion curves
m4_data <- mutate(m4_data,
  tor_diast = tf_rectify(tor_diastole),
  apm_diast = tf_rectify(apm_diastole),
  ppm_diast = tf_rectify(ppm_diastole)
)

m4_data %>% ggplot() +
  geom_spaghetti(aes(y = tor_diast)) +
  labs(title = "Rectified diastole curves") +
m4_data %>% ggplot() +
  geom_spaghetti(aes(y = apm_diast)) +
m4_data %>% ggplot() +
  geom_spaghetti(aes(y = ppm_diast))

Models

# data prep for pffr
m4_data$tor_d_mat <- m4_data$tor_diast |> as.matrix()
m4_data$apm_d <- with(m4_data, apm_diast - mean(apm_diast)) |> as.matrix()
m4_data$ppm_d <- with(m4_data, ppm_diast - mean(ppm_diast)) |> as.matrix()
# diastolic concurrent fof regression
m4 <- pffr(tor_d_mat ~ s(apm_d) + s(ppm_d) + Gender + Age + heart_rate,
  yind = seq(0, 1, length.out = 30),
  family = "gaulss", data = m4_data
)
summary(m4)
## 
## Family: gaulss 
## Link function: identity logb 
## 
## Formula:
## tor_d_mat ~ s(apm_d) + s(ppm_d) + Gender + Age + heart_rate
## 
## Constant coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)      0.783      0.729    1.07     0.28
## (Intercept).1    1.016      0.024   42.34   <2e-16
## 
## Smooth terms & functional coefficients:
##                      edf Ref.df Chi.sq p-value
## Intercept(yindex)  10.21  19.00   56.0 < 2e-16
## s(apm_d)           19.03  22.73  276.4 < 2e-16
## s(ppm_d)           13.48  17.23   60.2 2.4e-06
## Gender(yindex)      1.01   1.02    0.0       1
## Age(yindex)         3.18   3.54   49.9 < 2e-16
## heart_rate(yindex)  3.69   3.97   98.8 < 2e-16
## NA                 18.80  19.00 2264.5 < 2e-16
## 
## R-sq.(adj) =     Deviance explained =  100%
## -REML score = 5896.4  Scale est. = 1         n = 2250 (75 x 30)
pffr.check(m4, rep = 100)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 14 iterations.
## Gradient range [-0.0017,0.000295]
## (score 5896 & scale 1).
## Hessian positive definite, eigenvalue range [0.00169,6.2].
## Model rank =  125 / 125 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                             k'   edf k-index p-value
## s(yindex.vec)            19.00 10.21    0.86  <2e-16
## ti(apm_d,yindex.vec)     35.00 19.03    1.07       1
## ti(ppm_d,yindex.vec)     35.00 13.48    1.07       1
## s(yindex.vec):Gender      5.00  1.01    0.86  <2e-16
## s(yindex.vec):Age         5.00  3.18    0.86  <2e-16
## s(yindex.vec):heart_rate  5.00  3.69    0.86  <2e-16
## s.1(yindex.vec)          19.00 18.80    0.86  <2e-16
m4_gam <- m4
class(m4_gam) <- class(m4)[-1]

e <- getViz(m4_gam)

qq.gamViz(e, level = .9, CI = "quantile")

matrix(residuals(m4, reformat = FALSE), nrow = m4$pffr$nobs, byrow = T) |>
  cov() |>
  filled.contour()

plotly::plot_ly(
  z = ~ cov(matrix(residuals(m4, reformat = FALSE), nrow = m4$pffr$nobs, byrow = T)),
  type = "surface",
  colorscale = "Viridis", # Color scheme - you can change this
  showscale = TRUE # Show color scale legend
)
matrix(residuals(m4, reformat = FALSE), nrow = m4$pffr$nobs, byrow = T) |>
  cor() |>
  filled.contour()

plotly::plot_ly(
  z = ~ cor(matrix(residuals(m4, reformat = FALSE), nrow = m4$pffr$nobs, byrow = T)),
  type = "surface",
  colorscale = "Viridis", # Color scheme - you can change this
  showscale = TRUE # Show color scale legend
)
# changed color scheme here so we can easier distinguish between
# time-value regions with no significant association to TOR (white) and
# small-effect-but-still-statistically-significant regions (grey-ish hues)
plot(sm(e, 2)) +
  l_fitRaster(pTrans = function(.p) .p < 0.05) +
  scale_fill_gradient2(low = "blue", mid = "lightgrey", high = "red") +
  l_fitContour(bins = 20) + labs(y = "Diastole proportion", x = "APM", title = "") +
  guides(fill = guide_legend(title = "Effect on untwist"))
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

plot(m4, scheme = 1, select = 2, ticktype = "detailed")

plot(sm(e, 3)) +
  l_fitRaster(pTrans = function(.p) .p < 0.05) +
  scale_fill_gradient2(low = "blue", mid = "lightgrey", high = "red") +
  l_fitContour(bins = 20) + labs(y = "Diastole proportion", x = "PPM", title = "") +
  guides(fill = guide_legend(title = "Effect on untwist"))
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

plot(m4, scheme = 1, select = 3, ticktype = "detailed")

diastole_summary <- m4_data %>%
  pivot_longer(cols = 12:14, names_to = "Curve", values_to = "strain") %>%
  group_by(Curve) %>%
  summarise(
    strain_mean = mean(strain),
    strain_sd = sd(strain)
  )

diastole_summary %>% ggplot() +
  geom_spaghetti(aes(y = strain_mean, color = Curve)) +
  geom_errorband(aes(
    ymax = strain_mean + 2 * strain_sd / sqrt(nrow(m4_data)),
    ymin = strain_mean - 2 * strain_sd / sqrt(nrow(m4_data)),
    fill = Curve
  )) +
  labs(
    x = "Diastole proportion",
    y = "",
    caption = "mean+/-2 standard errors",
    legend = "Curve"
  ) +
  theme_bw()